Structure and Dynamics of a Phase-Separating Active Colloidal Fluid 
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We examine a minimal model for an active colloidal fluid in the form of self-propelled Brownian 
hard spheres that interact purely through excluded volume. Despite the absence of an aligning 
interaction, this system shows the signature behaviors of an active fluid, including anomalous number 
fluctuations and phase separation behavior. Using simulations and analytic modeling, we quantify 
the phase diagram and separation kinetics. The dense phase is a unique material that we call an 
active hexatic, which exhibits the structural signatures of a crystalline solid near the crystal- hexatic 
transition point, but the rheological and transport properties associated with a viscoelastic fluid. 
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Active fluids composed of self-propelled units occur in 
nature on many scales ranging from cytoskeletal fila- 
ments and bacterial suspensions to macroscopic entities 
such as insects, fish and birds [1]. From a physicist's 
perspective, these constitute a class of novel materials 
that exhibit strange and exciting phenomena such as 
dynamical self regulation [2], clustering [3], anomalous 
density fluctuations [4], and strange rheological behav- 
ior [5-7]. Motivated by these findings, recent experi- 
ments have focused on realizing active fluids in nonliving 
systems, using chemically propelled particles undergo- 
ing self-diffusophoresis [8-10], Janus particles undergoing 
thermophoresis [11, 12], as well as vibrated monolayers 
of granular particles [13-15]. 

In this letter we explore a minimal active fluid model: a 
two-dimensional system of self-propelled smooth spheres 
interacting by excluded volume alone. Unlike self- 
propelled rods [16-20], these particles cannot interchange 
angular momentum and thus lack a mutual alignment 
mechanism. Recent simulation and experimental stud- 
ies have shown that this system exhibits giant number 
fluctuations [21] and athermal phase separation [21, 22] 
that are characteristic of active fluids [4, 23, 24]. Here 
we employ extensive Brownian dynamics simulations to 
characterize the phase diagram of this system and its 
material properties. We find that the dense phase is a 
dynamic new form of material that we call an 'active hex- 
atic'. This material exhibits structural properties consis- 
tent with a 2D colloidal crystal near the crystal-hexatic 
transition point [25, 26], but is characterized by such 
anomalous features as superdiffusive transport at inter- 
mediate timescales and a highly heterogeneous and dy- 
namic stress distribution (see Fig. (1)). We also map the 
system's phase diagram over a large region of parameter 
space, develop an analytic model which captures its es- 
sential features, and draw an analogy between this ather- 
mal active phase separation and the behavior of equilib- 
rium systems with attractive interactions. The experi- 
mental accessibility, tunable phase behavior, and unique 
material properties of this system make it an attractive 
starting point for the design of new active materials. 

Model and Simulation Method: Our system consists of 
smooth spheres immersed in a Brownian solvent and con- 
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FIG. 1. (color online) A visual summary of the results of this 
work. Top left: Beyond critical density and activity levels 
the active colloidal fluid phase separates into dense active 
hexatic and dilute active fluid phases. The dense clusters 
coarsen until a single cluster remains (see SI in [27]). Top 
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right: The static structure factor S(k) = 

restricted to the interior of a large cluster. These signatures 
resemble those of a high temperature colloidal crystal near 
the crystal-hexatic phase transition. Bottom left: A heat 
map of the pressure inside the active hexatic. It is highly 
heterogeneous and dynamic, indicating a complex response 
to an external stress (see S2 in [27]). Bottom right: The 
mean square displacement of a tagged particle in the active 
hexatic. At intermediate time scales, it exhibits anomalous 
superdiffusive transport. 



fined to a plane. Each particle is self-propelled with 
a constant force, and interactions between particles re- 
sult from isotropic hard-core repulsion only. We include 
no mechanism for explicit alignment or transmission of 
torques between particles. 

The state of the system is represented by the positions 
and self-propulsion directions {vuOi\f =1 of all particles. 
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FIG. 2. (color online) Left: Phase densities as a function of 
Peclet number (Pe) for a range of overall <j). At low Pe the sys- 
tem is single-phase, while at increased Pe it phase-separates. 
The coexistence boundary is analogous to the binodal curve of 
an equilibrium fluid, with Pe acting as an attraction strength. 
Right: 'Giant' number fluctuations in high cj> and Pe systems. 
Plot shows the exponent 5 with the standard deviation of par- 
ticle counts in various-sized subsystems scaling as a ~ (n) s . 



FIG. 3. (color online) Left: Contour map of cluster frac- 
tion / c (Pe, (j)) measured from simulations. The dashed curve 
marks the approximate boundary of the phase-separated re- 
gion. Right: Cluster fraction as predicted by our analytic 
theory (Eq. 3), showing good agreement in the high Pe re- 
gion. The apparent disagreement for high 0, low Pe occurs 
because the identification of clusters becomes ambiguous near 
the random close packing density. 



Their evolution is governed by the coupled overdamped 
Langevin equations: 



Dp [F ex ({ ri }) + F pZ >i] 



(1) 
(2) 



Here F ex is an excluded- volume repulsive force given by 
the WCA potential V ex = 4k B T U^) 12 - (^) 6 j + k B T if 

r < 2^, and zero otherwise [28], with a the particle di- 
ameter. F p is the magnitude of the self-propulsion force, 
Vi = (cos#i, sin#i), and /3 = j^r* D and D r are trans- 
lational and rotational diffusion constants, which in the 
low- Reynolds-number regime are related by D r = 
The r] are Gaussian white noise variables with (rji(t)) = 
and (r h (t)r ]j (t f )) = 5 ij 5(t-t f ). 

We non-dimensionalized the equations of motion using 

2 

a and k B T as basic units of length and energy, and r = 
as the unit of time. Simulations employed the stochas- 
tic Runge-Kutta method [29] with maximum timestep 
2 x 10 -5 r. Simulations mapping the phase diagram were 
run with 15, 000 particles until time lOOr, while systems 
of 128,000 particles were used to explore kinetics and 
material properties. The simulation box was square with 
periodic boundaries, with side length chosen to achieve 
the desired density. The system is parametrized by two 
dimensionless values, the packing fraction <fi and the 
Peclet number Pe = v p ^, with the propulsion velocity 

from near-zero to 



and Pe from zero to 



v p = Df3F p . In this work, we varied 
the close-packing value </> CI 
150. 

Phase Separation and Large Number Fluctuations: We 
first show that our results are consistent with prior sim- 
ulations [21] and confirm that this system, despite the 
absence of aligning interactions, shows the signature be- 
haviors of an active fluid. In particular, the active spheres 
undergo nonequilibrium clustering (Fig. (1)) and exhibit 



giant number fluctuations (Fig. (2b)), similar to other 
model active systems [3, 19, 20, 30]. 

We next establish that this clustering is indeed ather- 
mal phase separation by measuring the density in each 
phase at different parameter values (Fig. (2a)). We iden- 
tify a critical point near (Pe = 50, <j) = 0.7) beyond which 
the system separates into two phases whose densities are 
independent of overall system density and are determined 
by the strength of activity alone. The result resembles 
the binodal curve of an equilibrium system of mutually 
attracting particles undergoing phase separation, with Pe 
playing the role of an attraction strength. The physi- 
cal mechanism underlying this phenomenon is a density- 
dependent suppression of the self-propulsion velocity (see 
Fig. S7 in [27]), leading to a negative correction to the 
macroscopic diffusion coefficient [21, 24, 31]. 
The Phase- Separated Steady State: To characterize the 
steady state, we measured the fraction of particles in the 
dense phase at time lOOr (Fig. (3)). In contrast with 
recent work [21] which reported a single critical density, 
we observe that this cluster fraction is a nontrivial func- 
tion of the system parameters / c (Pe, 0). To understand 
this relationship we develop a minimal model in which 
this function can be found analytically. Let us assume 
the steady state contains a macroscopic cluster which we 
take to be close-packed. Particles in the cluster are sta- 
tionary in space but their 0i continue to evolve diffusively. 
We treat the gas as homogeneous and isotropic, and as- 
sume that a particle colliding with the cluster surface is 
immediately absorbed. This model captures the physics 
of the density-dependent propulsion velocity [21, 24, 31] 
in a tractable manner (see Fig. S7 in [27]). 

Within this model, we can write the rate of absorp- 
tion of particles of orientation from the gas phase as 
&in(0) = ^~Pg%> cos #? where p g is the gas number den- 
sity. Integrating yields the total incoming flux per unit 
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FIG. 4. (color online) Left: Defect structures in a large clus- 
ter. Regions of high crystalline order (white) coexist with 
isolated and linear defects (dark). The color of each particle 
indicates its \qe\. Inset shows pairs of 5/7 defects (red/blue). 
Right: Correlation function (qg (r)qe(r')) for clusters at vari- 
ous Peclet numbers, showing a transition from exponential to 
power-law decay as activity is increased. At high Pe, decay is 
slower than |r — y'\~^ (red dashed line). 

length: h in = ^2.. To estimate the rate of evaporation, 
note that a particle on the cluster surface will remain 
there so long as its self-propulsion direction remains "be- 
low the horizon", i.e., n • i> < 0, where n is normal to 
the surface. When its direction moves above the horizon, 
it immediately escapes and joins the gas. This rate can 
be calculated by solving the diffusion equation in angu- 
lar space with absorbing boundaries (for clusters large 
enough to treat the interface as flat, at ±|) and ini- 
tial condition given by the distribution of incident par- 
ticles: d t P(0,t) = D T d$P(0,t), with P(±f,t) = and 
P(#, 0) = \ cos 6. Further, the departure of a surface 
particle creates a hole through which subsurface parti- 
cles (whose t>i may point outwards) can escape. With k 
we denote the average total number of particles lost per 
escape event, which we treat as a fitting parameter. The 
total outgoing rate is then fc out = ^ JL - 

Equating h in and fc ou t yields a steady-state condition 
for the gas density: p g = . p g can be eliminated 
in favor of / c , yielding (in terms of our dimensionless 
parameters): 

= ^e-y* (3) 

This function is plotted in Fig. (3) with k = 4.5, in good 
accord with our simulation results. Further, the condi- 
tion f c = allows us to deduce a criterion for the onset 
of clustering. Restoring dimensional quantities, this con- 
dition gives 4>av p ~ D r . Note that (j)crv p is a collision 
frequency; thus the system begins to cluster at param- 
eters for which the collision time becomes shorter than 
the rotational diffusion time. 

Structure of the Dense Phase: Since the system is com- 
posed of monodisperse spheres, the dense phase exhibits 
crystallization [32]. As shown in Fig. (1) the static struc- 



ture factor of the cluster interior shows well-developed 
sixfold symmetry when activity is high. Further, the ra- 
dial distribution function shows clear peaks at the sites 
of a hexagonal lattice (see Fig. S6 in [27]). The num- 
ber and sharpness of the peaks increase with Pe. We 
also measured the bond-orientational order parameter 

26 W = WUT\ ^ieAT(i) ei69ij i where -^W runs over the 
neighbors of particle i (denned as being closer than a 
threshold distance), and Oij is the angle between the 
i-j bond and an arbitrary axis (Fig. (4)). We find a 
structure characterized by large regions of high order 
with embedded defects that are predominantly 5-7 pairs 
(Fig. (4a) inset and S4 in [27]). Next, we measured the 
correlation function (q% (r)qe(r f )) as a function of \r — r'\ 
(Fig. 4), finding a power-law decay for systems of high 
activity with exponents < 1/4, which puts the system 
above the hexatic-liquid coexistence region as predicted 
by KTHNY theory [33]. Based on this evidence, the sys- 
tem is near the coexistence regime of a crystal and a 
hexatic. This active material is unique in that it is held 
together by active forces alone and the consequent arrest 
of motion due to frustration. In this sense it is similar to 
amorphous materials such as granular packs as reflected 
by the highly heterogeneous stress distribution (Fig. (1)) 
[34]. 

Dynamics in the Dense Phase: Within the active hex- 
atic material, self-propulsion forces continuously evolve 
by rotational diffusion, breaking local force balance and 
leading to defect formation and migration (see S4 in [27]). 
A compelling way to view the motion produced by this 
athermal process is a simulated FRAP experiment [35], 
in which particles within a contiguous region are tagged, 
making subsequent mingling of tagged and untagged par- 
ticles visible (see S3 in [27]). To quantify this behavior, 
we measured the mean square displacement (MSD) of 
particles in the cluster interior. As shown in Fig. (1), we 
observe subdiffusive motion on short timescales, followed 
by a superdiffusive regime, returning to diffusive motion 
on long timescales. The exponents of the subdiffusive and 
superdiffusive motion (| and |, respectively) are well- 
conserved across a wide range of propulsion strengths. 

The observed dynamics within the active hexatic can 
be viewed in two ways. First, note that an isolated self- 
propelled particle will exhibit diffusive, ballistic and dif- 
fusive behavior on time scales t<^-,^-<t<-^- and 

t > jy- respectively (see Fig. S5 in [27]). These dynamical 
regimes are modified by the active hexatic environment; 
in particular, the ballistic regime is modulated by "stick- 
ing" events as the particle is localized in crystal domains, 
resulting in the observed Levy flight like behavior [36, 37]. 
Second, treating the observed MSD behavior as a window 
into the material's rheology [38] enables us to use the re- 
lationship between the shear modulus and particle MSD 
in Laplace space, G(s) = 7rscr (A^(5)) •> ^° see ^ na ^ a ) a ^ 
short time scales (high frequency), the effective viscosity 
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FIG. 5. (color online) Examples of phase separation kinet- 
ics. Left: A system with Pe = 100, <j) — 0.45 in which a 
delayed nucleation event leads quickly to steady-state. For 
shallowly- quenched systems, the nucleation time can be long 
enough that artificial seeding is needed to make nucleation 
computationally accessible. Right: A system with Pe = 80, 
cj) — 0.6 where spinodal decomposition leads to a coarsening 
regime which slowly evolves towards steady-state (see SI in 
[27]). Inset shows mean cluster size scaling approximately as 
ti. (see S8 in [27]). 

scales as cj -1 / 2 and hence the material is shear thinning, 
b) at intermediate frequencies, the viscosity scales as a; 1 / 2 
and hence is shear thickening and c) on the longest time 
scales the material is Newtonian in its response to shear. 
Thus, the rotational diffusion time, which sets the per- 
sistence of motion in the gas phase, also demarcates the 
transition from complex to Newtonian rheology in the 
active hexatic material. 

Kinetics of Phase Separation: Despite the at hernial 
origins of phase separation in this system, simula- 
tions quenched from a homogeneous state to parame- 
ters where clustering occurs move through familiar nu- 
cleation, growth, and coarsening stages (Fig. (5)). How- 
ever, in the coarsening regime we find surprisingly that 
the mean cluster size scales as with a correspond- 
ing length scale C(t) ~ t* (Fig. (5) inset, also see S8 
in [27]). This differs from the standard 2D coarsening 
exponents, but matches recent simulation results for the 
Vicsek model and related active systems [39] . Our results 
should be viewed as preliminary due to the limited range 
of our data, but nevertheless this unexpected similar- 
ity between the coarsening of point-particles with polar 
alignment and that of hard spheres with no alignment 
suggests a deep relationship between these very different 
types of systems. Future work is needed to uncover the 
origins of these scaling exponents and their implications 
for universality in active fluids. 

Discussion: A fluid of self-propelled colloidal spheres 
exhibits the athermal phase separation that is intrinsic 
to active fluids and is a primary mechanism leading to 
emergent structures in diverse systems [2, 24]. We have 
shown that the physics underlying this phase behavior 
can be understood in terms of microscopic parameters. 
From a practical perspective, our simulations show that 



the active hexatic dense phase exhibits a combination of 
structural and transport properties not achievable in a 
traditional passive material. Further development of ex- 
perimental realizations of this system (e.g. Ref. [22]) 
will advance the development of materials whose phase 
behavior, rheology, and transport properties can be pre- 
cisely controlled by activity level. 
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